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The manuscript establishes a series expansion of the core integral that relates changes in longitude 
and latitude along the geodetic line in oblate elliptical coordinates, and of a companion integral which 
is the path length along this line as a function of latitude. The expansion is a power series in the 
scaled (constant) altitude of the trajectory over the surface of the ellipsoid. Each term of this series 
is reduced to sums over inverse trigonometric functions, square roots and Elliptic Integrals. The 
aim is to avoid purely numerical means of integration. 
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I. SCOPE 
OA A. Geodetic Coordinates 



.^ , An ellipsoid is a reference surface fixed by an equatorial radius pe and a polar radius pe ■ In many applications the 

u ■ 

^': pl^pl{l~e% (1) 

-(— » ■ 

■ is the principal reduced parameter. The three-dimensional ellipsoidal coorcfinatcs, altitude h and the angles of 

longitude A and latitude <j>, are basically defined with the aid of a straight plumb line along the shortest distance 
between a general point and its "foot" point on the surface. The relation between the Cartesian geocentric coordinates 
{x,y,z) and the curvilinear {h,X,(p) is '6, 8 10, 12 14], 
> 

O' ( ^\ ( [N {(l)) + h] COS (j) COS X 

0\\ \y\ = \ [N{4') + h]cos(j)siTi\ |, (2) 

r^ ; \z ) \[N{(j)){l-e^) + h\sm.(j) 

l/~) where 

o : N{<j>) ^ -; (3) 

— ' V 1 — e^ sm (j) 

is the distance between the foot point and the polar axis measured along the straight extension of the plumb line. 
Jv>( , The projection t onto the polar axis, 

$^ 

C^ ■ T = sin0, (4) 

will be useful to substitute trigonometric functions by rational functions. 

B. Geodetic Lines 

A geodetic line is the line of shortest Euclidean distance between two points within a surface of constant height 
h. This balances differential changes in the trajectory 4){X) with the two principal curvatures at each point; in 
consequence, the meridional radius of curvature 
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often appears to condense the notation. 

The solution of the differential equations of the geodetic lines crystallizes in the integral [llj 



I{t)= f dr 



cih + M) 



{N + hf{l-T-^)Jl-T 



= AA, 



(6) 



(N+hy 



which relates a difference in latitude — the limits of the integral — to a difference in longitude — the right hand side. 
The obliquity parameter c picks an individual geodetic line out of the bundle of all lines that cross a general point. 
Considering the r at which the discriminant of the root in the denominator of (jH]) is zero shows that c is also the 
distance to the polar axis at the point highest above (or below) the equatorial plane .llj . 

The single interest of this manuscript is in demonstrating a semi-numerical approach to evaluation of this integral. 
The constant of integration is tacitly fixed to imply the lower limit t — 4) — Qm the integral, because such a reference 
to the nodal line leads to well defined branch cuts of all square roots involved. The strategy is to expand the integrand 
into a power series of /i/pe, and to exchange the order of integration and summation, which defines a family of integrals 
with two additional parameters reminiscent of the order in the expansion. Each of these is reduced to the level of 
multiple — but finite — sums over Elliptic Integrals, assuming that these are accessible through a numerical library 

In overview, one way of computation of © is addressed in Section |lTl Auxiliary integrals fall into two classes, one 
reducible to roots and inverse trigonometric functions, the other to elliptic integrals. The distance along the geodetic 
line defines another integral which is treated in the same spirit in Section IIIII Its power series yields a family of 
integrals which can be recast for efficient reuse of the functionality build in Section |lll Finding the inverse of / with 
respect to the parameter c is closely related to the inverse problem of geodesy and shortly addressed in Section |lVl 



II. LONGITUDE-LATITUDE COUPLING INTEGRAL 



A. Taylor Expansion in Powers of Altitude 



Expansion of the auxiliary N and M and lifting of some square roots provides a long write-up of ([5]), 



dr- 



A(i„eV2)3/2 + l_e2 



(1 + AVl-e2T2)(l - T2)Vl-e2T2./(l + ±-^l^e^T^y{l _ ^2) - (c/pe)2(l - c'^t'^) 



(7) 



which we intend to calculate. The altitude h and parameter c appear only scaled with pe, so introducing a function 
of two dimensionless variables h and c. 



Iaih,c) = dr 



{l-e^T^y 



(1 + Wl-e2r2)(l - r2)y^(I + hVl ~ e^T^)^{l - t^) - c\l - e^r"^) 



shows the composition 



Pe 



-/i(-, -) + (1 - e2)/_i/2(-, -) 

Pe Pe Pe Pe Pe 



The structure of the integrand is dominated by the two variables 



T=1-t'^\ £' = I-eV^ 



The power series of (|8]) becomes 



Ia{h, c) = ^{-hy ^ K,,fe / dr 

c- — n 7 n '^ 



j^a+s/2 



ifc-l 



(r-c2£;)fc+i/2- 



(8) 



(9) 



(10) 



(11) 



s=0 fe=0 

The coefficients k emerge from the product of the geometric series of 1/(1 + hE^^^) by the binomial expansion of 



l/y^{l + hE^/^yT-c'^E in 



Ks.k = 4' 



l-k 



E 

l=k 



(-i/2)' = (-i)M '/') + (-i)^2^'=— 1( '/^)Hl-:--'r'Ho)- (12) 
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TABLE I: Table of the rational values of Ks,fe by equation (|12|l . Each column attains a constant value for rows s > 2k. 
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FIG. 1: Values of /_i/2,fc with a lower limit of r = for three different c and three different k with e = 0.08182. Small c 
indicates near-polar routes where the branch cut with Ip,k — >■ oo is at larger r. 



The term with the Jacob! Polynomial P is to be interpreted as zero if fc < s/2. (fTT|) states the problem in terms of 
integrals 



//3^fc = / dr 



I' 



rpk-l^[3 



[T - c2^)'=+l/2 



= / dT- 



(l-r2)'=-i(l 



= 2^2)/3 



g2^2^1*:+l/2■ 



(13) 



for P = a + s/2 = —1/2, 0, 1/2, 1, 3/2, . . . and fc = 0, 1, 2, 3, . . .. For small eccentricities, the Ip^k are nearly independent 
of j3 because E is close to unity, so it suffices to illustrate the values for zero lower limit and variable upper limit r 
for one value of /3 in Figure [TJ The main disadvantage of the method is that the eventual constancy of K^^k down the 
columns of Table |I] in conjunction with the alternating sign of {—h-Y in (jlip induces oscillatory behavior (cancellation 
effects) of the series if h is not small. 



For integrals classifications, (TTSl) is phrased as 



(l-T2)fc-l(a2-T2)/5 



where two parameters a and 6, 



-^/3.fe (1 _ c2e2)fe+l/2 Jg "^^ (62 _ 7-2)fe+l/2 



a=l/e2>6 = (l-c2)/(l-c2e2) 



(14) 



(15) 



decide on the branches of the Elliptic Integrals. The exponent /3 either leads to reduction to elementary functions if 
it is an integer (Section III Bp . or to Elliptic Integrals if it is half- integer (Section III C[) . 



B. Cases of Elementary Functions 

If /3 = 0, 1, 2, .. and fc > 0, we substitute x = r^ in (ITi|) . 

^'^ ~ 2(1 - c2e2)'=+i/2 jg '^ ^(52 _ a;)fc+i/2 



Jl3,k{x) ^ dx 



k-ifn"^ -T.^/3 



{l-xY-\a^ ~x) 



equivalent to the computation of 



Vx(62 _ 2;)fe+l/2 

Partial integration of this integral generates the recurrence 

2(fc- 1)62 + 2/2 + 1 



1 + 2/3 + 



62-1 



Jp^k - 2Vx ^^2 _ 2;)fc+l/2 + 2^j— y J/3,fc-l 

62 

+2a2/3J^_i^fc + (2fc + 1)^^— - J^ ;,+!. 



This allows to build the entire table of Jpj. from a list of J^.o, >//3,i and Jg.fe. 
The substitution 

1 ,2 1 
z = — , X = 6 

r ~ X z 

and binomial expansion establish 

[(l-62)z + l]'=-i[(a2-62)^ + l]/3 



J, 



/3,fe 



dz- 



V62z - 1 






^s+m- 



V62z-1 



fc>l, /3>0. 



• The special case s + m — /3 = is covered by 



dz 



\/cz — 1 c 
• The cases s + m - /3 < are handled by f?,, 2.245.2] 



-\/cz — 1. 



dz tr(i + l/2) 



z*+iV^7^T r(i + i) 



Vcz - 1 ^ 



r(^ + i) i_ 

r(3/2 + (cz)'+ 



-j- + 2^^ arctan ycz^-T 



(16) 



(17) 



(18) 

(19) 

(20) 
(21) 

(22) 
t > 0. (23) 



The sum evaluates to zero if the upper limit is smaller than the lower limit. 
• The cases s + m — /3 >0 are solved by [7}, 2.222] 



Vcz — 1 



dz 



2y/cz - 1 ^ /A (cz - 1)' 
1=0 ^ ' 



t > 0. 



(24) 



This has effectively written (|2T]) as triple sums. [In numerical practice, the integral in (1211) is placed into a look-up 
table for the j3 + k different values of the exponent s + m — /3.] The case (3 — appears as a double sum, but is recast 
into a single sum by resummation of the l-sum in (|24|) and the s-sum in (f2T|) : 



fc-i 



1 /7^^ v-^A-A. ,9^;(&^^-l)' 



1=0 



I 



21 + 1 



2^^^2FiQ>i-fc;|(&^-i)(&'--i) 



52/c^ 



-.2F1 {l,l~k 





(25) 


3 (l-62)x\ 
'2' b^-x / ■ 


(26) 



The remaining part of Section IIIBI considers the case fc = which is not available from ([21 



=2/3 



^^,0 — 



2%/r 



c e 



dT 



2 _ ^2^/? 



(a2 - r2) 



(1-t2)VP" 



o2/3 



2Vl - c2e2 



i:0("^-i)'-'/ 



dT 



(l-r2) 



2\i-l 



V6^^^ 



„2/3 



2Vl - c2e2 



(a^ - 1)^ 



1 62(1+t2)2t2 , ^^/3 



2vr^^ 



E 



62(1 -r2) -^ VV («^-l) 



dT 



(l-r2) 



2\;-i 



\/62-T- 



fi2 _ ^2 



(27) 



[The split value at / = is derived substituting y = 1/(1 — r^), then using [7, 2.261]. The constant of integration is 

chosen to realize the limit /o,o — > 0.] 
This is a superposition of 



dT 



(l-r2) 



2\l-l 



;-i 



V62 - r2 -^-^ \ t 



^ ( \. ^ ) (1 - 62)'-!-* / dT(62 - T^^^- 



f dTib" - t2)'-1/2 - ^ [' ^ ^] (1 - h')'-^-'B,, 



(28) 



where we have defined 



Starting from 



B, = JdTib^-Ty-^/^; i = 0,1,2,., 



(29) 



the recurrence [^, 2.260.2] 



Bq = arcsin — , 




B.^l-Xb^-T^y-^f^+(i-l-)m^. 



(30) 



(31) 



allows to calculate (p8)) and eventually vG7\ 



C. Elliptic case 

This section looks at ([H]) for the cases /? = — -i, -1, |, | Let 



J^,fc(r) = 2 / dT 



(l_^2)fc-l(^2_^2)/3 
(-52 _ ^2-)fc+l/2 ■ 



(32) 



The factor 2 in the definition is chosen to maintain the "alignment" Jp^k{T) = Ji3,k{T'^)- Partial integration offers the 
recurrence 



1 + 2^ + 



2{k - l)b^ + 2k + 1 
62-1 






(62 - r2)fe+i/2 
+2a2/3 J^_i,fe + (2fc + 1) 



62-1 



62-1 



J. 



I3,k+1, 



(33) 



which is the same as p8)) . 

For fc > 1, binomial expansion of the numerator in (1321) proposes 



fc_l/3+l/2 



jp,k - 2 51 XI 



fc-1 



(l-62)^(^ + ^/^)(a2_52)'" 



dr 



These elhptic integrals with integer exponents u = 1 + s — (/3 + 1/2 — to) > are 3, 219.06] 

1 1 



dr 



^{a^^T^W~T^){b^ _ ^2Y ab^^ 



D2vi 



where ^, 313.05] 0, 3.158.15] 



Do = F(e,6/a) 
D2 = F{^,b/a) 



fe2 






i?(e,6/«) 



(34) 

(35) 

(36) 

(37) 



with sin^ = rjb. [In this equation and where used with two arguments, (|45p . (|56p and the last equation in the 
appendix, E is the incomplete Elliptic Integral of the second kind, elsewhere the shorthand ([TU|.] The recurrence 
extending these two values both sides to larger or negative v is 0, 313.05] 



{2v + l)(a2 - b')D2-„^2 = (2w - l)6'i?2.»2 + 2i;(a2 - 2b')D2-, + aV^^^^^^- 

cos^^ ^ 



(38) 



Similar to the exception in Section IIIBl the case fc = is not covered by the expansion above and established 
individually: 



J/3,0 = 2 / dr 



2_^2^/^ 



{a'-T^) 



(l-r2)^/52 



2 / dr- 



(a2 - t2)'3+1/2 



(1 - r2)^(a2-T2)(62_r2) 

^0 V "^ ^ 7 V(a2-r2)(62_^2) 

The term m = is an Elliptic Integral of the third kind d 219.02] 0, 3.157.7]: 

dT , ^ = = -^(C,6^6/a). 

(1-T2)^(a2-T2)(fe2„r2) a ^^' 



(39) 
(40) 

(41) 
(42) 



Because the case fc = s = with /3 = —1/2 is the only contribution to pT|) and ([5]) if the geodetic line is on the 
surface of the ellipsoid (h — Q\ this is the only value relevant to the integral (O for this "classic" case. 
The terms to > 1 are delegated to [3, 219.05] by binomial expansion of the numerator. 



dT- 



(l-r2) 



2\m-\ 



^{a^-T^){b^-T^) 



E 

1=0 



I 



2\l 



TO - 1\ (-62) 



12; • 



starting from [3, 310.02,310.05] 



Ao = F{^,b/a); 
b^A2 - a^[F{^,b/a)~E{^,b/a)] 



(43) 



(44) 
(45) 



more values of b^^A2i may be generated from the recurrence [3|, 310.05] 



{21 + l)b''+^A2i+2 = a^b-^-T-^^a'^-T^ + 2l{a' + b')b''A2i + (1 - 2l)a'b'' A21-2; I > -1- 



(46) 



III. LINE DISTANCE INTEGRAL 

A. Reduction to the Angular Coupling Integral 

The formula for the distance s along the geodetic line is given by [11[ 



(l-e2T2)./l-T2- 



w+w 



Pe I dr'- 2^ ^"^ , '^\ ' (48) 



= P: 
calling a group of integrals 



h ,h c , , h ^^ , h c ^ , „, , h c ^ h , o, , h c 



Soi-, -) + i-rSi/2i-, -) + (!- e^)5-3/2(-, -) + -(!- e')S^ii- 

Pe Pe Pe Pe Pe Pe Pe Pe Pe Pe Pe 



(49) 



So.{h,c)= f dr , , ^" (50) 

with two dimensionless scaled parameters and one characteristic exponent a. Binomial expansion of the square root 
provides a power series in h, 

s=0 k=\s/2] \ / \ / J \ I 

The integral is similar to ([13]); the difference is a factor T plus the request from (|49)) to evaluate the cases a — —3/2 
and a = — 1. The factor T is distributed with the aid of 

T = ^-^ + ^,E. (52) 

Recalling p3|) . this maps ([ST]) on the integrals p^ 

'^'^ (T - 2fJ\k+l/2 ^ (-*- ^ 0,)Ia+s/2,k + a-^a+l+s/2,fc- (53) 

B. Special Values 

To carry out the right hand side of ([55]) . the only aspect not yet covered by Section |ll] is to implement I-3/2,k ^^nd 
I-i^k- Furthermore, k is only required for the restricted range of s seen in the summation (|5ip . which reduces the 
"new" cases further to 

• -^-3/2.0 from (s, fc,a) = (0,0,-3/2), 

• /-i,o from {s,k,a) — (0,0,-1), 

• and /-i,i from (s. fc, a) — (1, 1, —3/2), 

because otherwise the first index of Ij3^k is /? > —1/2, already treated in Section |lTl Turning to 

^-3/2,0 = (i_^2g2)l/2 J ^^ (1 - r2)(a2 - r2)3/2(62 _ ^2)1/2 (^4) 

proposes partial fraction decomposition 

dr- „, , „ „,„,„, „ „, , ,„ = -^ / dr- 



(1-T2)(a2_ ^2)3/2(52 „ ^2)1/2 a2 - 1 7 (1 - r2 ) ^ (a^ - r2) (62 _ ^2) 



^ ^dr , ^ =, (55) 

(a2-r2)^(a2-T2)(62_r2' 



and the two integrals on the right hand side are known [3, 219.02,219.07,315.02] 

1 f„.. .o . , , 1 



-3/2,0 



e(a2 - l)Vl - c2e2 



m,b^b/a) 



b^ 



E{C,b/a)~- 



T Vb^ - r2 



a y/a'^ - r^ 



Demonstrated in (j49| . this is the only 5*0 value required on the surface of the ellipsoid — where h — 0. 
The second remaining case is 



/_ 



1,0 



dr- 



1 



VI - C2e2 J (l-T2)(a2-T2)V62 



which splits into two partial fractions — equivalent to ([55)) — with known integrals [3, 2.284] 

1 



/-i,o 



VI - c2e2(a2 - 1) 
The third remaining case is 



VT^F 



: arctan r 



l-fe2 

62 -t2 



1 /r /a2-62 

— , arctan [ —\ -r^ ^ 

aVa2 - fe2 I a V fe2 - r2 



/-i.i = 



dr 



(l-c2e2)3/2 7 (a2_^2)(^2_ ^2)3/2' 



with partial fractions 



dr 



dr 



(a2-T2)(62_r2)3/2 b^-a^J (^2 _ ^2y^2 _ ^2 b^-a^J (fe2_^2)3/ 



dr 



(56) 



(57) 



(58) 



(59) 



(60) 



The first integral on the right hand side is the same as met while calculating /-i,o- The second is also known [3, 
2.271.5], to yield 



I-hi = 



(l-c2e2)3/2(a2_52) 



b^V^ 



Vo 



62 



arctan 



T a? — b"^ 



b^-T^ 



(61) 



IV. INVERSE FUNCTION 



In sections HIl and Hill the value of the parameter c was regarded as known. The "inverse" problem of geodesy, on the 
other hand, is finding c assuming the value of the integral, i.e., AA, and of its limits are given. If Newton methods are 
employed to this problem, they also call for computation of the power series of dl/dc, which is addressed as follows: 

The derivative of (O is 



dj = -/ + - 

Pe Pe 



^dMh,c) + -(1 - e^)dj^i/2{h,c) 

Pe Pe 



by the product rule. The derivative of (|lip is given by the chain rule: 

oo s „ 



j^l+a+s/2 

[T - c^Ef+'^l'- 



:T 



fc-1 



(62) 



(63) 



The exponent of T in this integrand is deficient by 1 compared to the definition of Ip^k in (|13p — whereas it is abundant 
by 1 in (j5ip . Still, new functionality is not required, because 



(T - c2^)'=+3/2 



dr 



j^l+a+s/2 



T{T - c?E) {T - d^Ef+^l^ 



rpk 



(64) 



= / dT- 



c?E 



T-c^E T 



^l+a+s/2 I 1 

rpk _ _T _T 



(65) 



converts to integrals already discussed in Section |TT1 



V. SUMMARY 

Given the altitude, a directional parameter and a starting position, the task of finding the trajectory of the geodetic 
line in 3-dimensional geodetic coordinates turns into the evaluation of an integral which emerges from the solution 
of a differential equation which couples longitude and latitude. The coefficients of a systematic expansion of this 
integral in a power series of the altitude (scaled by the equatorial radius) have been reduced to multiple sums over 
elementary functions and incomplete Elliptic Integrals; the geodetic line on the surface of the ellipsoid is a specific 
and the simplest case. Auxiliary integrals developed along this path recur if the same expansion strategy is applied 
to other integrals related to the geodetic line. 

Appendix A: Table Errata 

Errata to the 1981 edition of the Tables of Sums, Products and Integrals [7] relevant to this work are: 

• 2.284: Preserve the sign of c on the right hand side by moving it out of the square root: 

Ax + B , A 2Bc-Ab 

zdx = — /i H . ^ J2- 



ip + R)y/R c c^p[b^~4.[a + p)c] 



• 2.245.2: Alternate signs of both b on the right hand side: 

"2*^-1(71- l)(n-2)---(n-fc)A'^ t 



^^ (2n-2m-3)(2n-2m-5)---(2n-2m-2/c + l)(-6)'=-i 1 , 

k=2 

(271 - 2m - 3)(2n - 2to - 5) • • • (-2m + 3)(-2m + l)(-6)"-i f z'^dx 



2"-i(n-l)!A" J t^i 



• 



3.158.15: Remove a & in a numerator of the right hand side: 

dx 1 „. ^ 1 



yJia'^-x'^W-x'^Y ab^ '" ' b^a^-b'^) 



X { aE{r], t) - u\l ^ ^ \ , [a > 5 > u > 0]. 
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